The broadcasting of multivariate RVs is very weird already at the Aesara level:

aesara.tensor.random.basic.multivariate_normal(
    mean=np.ones((3, 9)), cov=np.eye(9), size=(3,)
).eval().shape
>>> expected: (3, 9)
>>> actual: (3, 3, 9)

The above parametrization apparently gave it a 2-dimensional support? That's a problem, because the ndim_supp attribute belongs to the Op and is still at 1.
The other interpretation - keeping support fixed at 1 - is that it decided to replicate twice? Also doesn't make a lot of sense to me.

I noticed this when I tried to refactor the Distribution.dist() code such that size is eagerly passed instead of creating rv_native with size=None.
I attempted this refactor because it seemed more elegant, but in the resize-based implementation seems to be way more robust.

This can be rebased onto #4693 and, instead, re-introduce a non-breaking version of the shape keyword.

Closing in favor of #4696.